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NONLINEAR INTERACTION MODELING 
IN PHOTONIC CRYSTALS 


Dana Georgeta POPESCU,’ Paul STERIAN? 


Abstract. This paper presents several:methods for photonic crystals (PCs) characterization 
which are subsequently applied to realistic systems using specialized software. These methods 
are Finite Element Method, Plane Wave Method, Finite Difference Time Domain Method, 
Multipole and Bloch Wave Methods. Based on a full nonlinear system of equations, a 
nonlinear SHG (Second Harmonic Generation) problem in one-dimensional photonic crystals 
will be analyzed using Finite Element Method coupled with Fixed Point Iterations. The 
diagrams resulted after parameter variations are plotted using FlexPDE Professional 
software. As a result we obtained two maximum intensities of the second harmonic wave 
within each high index layer, contrasting to the fundamental wave peak. Also a representative 
band structure using the Plane Wave Method is given. 
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1. Introduction 


Photonic crystals are periodic dielectric structures which have a band gap that 
forbids the propagation of radiation within a certain frequency range. This 
property allows an enhanced control of light effects, otherwise difficult to control 
with conventional optics. 


Photonic crystals are described by Maxwell’s equations: 


BW (1.1) 
V-B=0 (1.2) 
vxE=-% (1.3) 

Vx B= fl #6 (1.4) 
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where E is the electric field strength, B is the magnetic flux density, p is the 
charge density, ¢, and sz, are respectively the vacuum dielectric constant and 


permeability and J is the current density. 


We can classify this materials with *photonic band gaps” (PBGs) in one-, two- 
and three dimensional [1], as illustrated in Figure 1.1. 


1-D 2-D 


1D periodic 2D periodic 3D periodic 


Fig. 1.1. Examples of one-, two, three-dimensional PCs. 
Different colors represent materials with different dielectric constant. 


The main goal of this paper is to present some theoretical models: Finite Element 
Method, Finite Difference Time Domain Method, Plane Wave Method, Multipole 
Method and Bloch Wave Method and then verify some of them for realistic cases. 
These techniques are very useful in understanding the special properties of the 
photonic crystals. In this manner we can predict and explain their behavior 
without resorting to brute force calculation. At the end, we highlight the utility of 
FlexPDE Professional and Matlab software in 1D and 2D photonic crystal 
computer simulations. 


2. Methods 


2.1. Finite element method in photonic crystals study (FEM) 


Using finite element method we can effectively investigate the modes allowed by 
photonic crystals when imposing periodic boundary conditions to a unit cell of a 
given material. Such PBG (photonic band gap) materials, restricted at the two- 
dimensional periodicity, have been successfully fabricated. We can obtain a 
reduction in complexity (a reduction in processor and memory requirement) for 
the in-plane propagation modes of photonic crystals using 2D Lagrangian finite 
elements. The FEM is efficient due to scattered eigensystem matrices and to the 
discontinuous dielectric constant which is handled in the real space. As an 
important application of the FEM we can consider the computation of the band 
structure for many common PC structures. With this numerical FEM technique we 
can easily and exactly represent the sharp discontinuities in dielectric constant and 
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solve boundary value problems for many PBG materials. The FEM is also used in 
meshing the dielectric function profile of the photonic crystals and by solving 
eigenvalue equation with proper periodic boundary conditions following the 
Bloch theorem we can calculate the in-plane band structures and in-plane photon 
density of states (PDOS) of the two-dimensional photonic crystals for the TE and 
TM modes. 


Referring to PC structure the FEM involves. breaking down the problem domain 
into many small elements of simple shape. The equations used are derived from 
Maxwell’s equations that approximate the electromagnetic field over an element. 
An interpolation function coupled with its coefficients express this approximation 
of the field. Making a compromise between the quality of approximation and the 
number of degrees of freedom (the number of function coefficients) we can 
choose an interpolation function. A good approximation of the true solution might 
be given by a higher order interpolation function, but this will require a large 
number of function coefficients thus increasing the computational and storage 
costs. The function coefficients must be computed for every element after 
choosing the interpolation function. These coefficients are then stored as 
elemental matrices, which may be subsequently assembled into global matrices by 
mapping local to global interpolatory node numbers. So an eigensystem of 
scattered matrices is formed by these global matrices. Using a subspace iterative 
technique we can solve this system of equations. 


The resulting eigenvalues are the frequencies of the allowed modes with their 
corresponding eigenvectors representing the field strength at the nodes [2-3]. 


A highly interesting case is prompted by introducing a defect into an otherwise 
highly ordered PC structure. The defect could be an alteration of a rods diameter, a 
variation in the dielectric constant of a rod or complete removal of a rod or rods [4]. 
The FEM output can certify the effect of these defects, used in controlling the 
frequency and range of photonic band-gaps. 


The modes having a frequency within the range of the band gap can propagate 
when we introduce a defect into the crystal structure. The defect induced mode is 
effectively “pinned” by the defect [1] and it cannot propagate through the crystal 
as it has a frequency within the forbidden range (the light is localized and cannot 
escape’). 


Generally speaking, the finite element method is a numerical technique used in 
finding approximate solutions of integral equations, or of partial differential 
equations (PDE). The solution approached is based on transforming the PDE into 
an approximate system of ordinary differential equations, or on completely 
eliminating the differential equation. This system is then integrated numerically 
using standard techniques such as Euler's method. 
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In solving PDE, the main challenge consist is finding an equation that 
approximates the original one. Great care should be paid in order to avoid the 
accumulation of errors from intermediate calculations and from the input which 
could make the resulting output meaningless. 


We can use FEM for solving PDE over complicated domains, when the desired 
precision varies over the entire domain, or when the domain changes. 


FEM is used in aeronautical, automotive industries, biomechanical, in design and 
development of some products, in minimizing materials, weight and costs. It 
allows the detailed visualization of the regions where structures twist or bend and 
additionally indicates the distribution of displacements and stresses. Also one can 
construct an entire, refined and optimized design before it is manufactured. FEM 
design has improved the methodology of the design process and the standard of 
engineering in many industrial applications [5-9]. 


FEM benefits include enhanced design, increased accuracy and_ better 
understanding of design parameters and other cycles in production. 


To illustrate the FEM we considered a nonlinear photonic crystal which offers 
unique and fundamental methods of enhancing various nonlinear optical 
processes. In the third part of this paper we used nonlinear Helmholtz equations, 
which are based on Maxwell's equations, and the corresponding equations of 
boundary conditions to emphasize the second harmonic effect. The nonlinear 
problem revealed has a unique solution only if the contribution of the nonlinear 
susceptibility tensor is not too large. We show that second harmonic generation 
efficiency can be improved by creating defects in the periodic structure. 


2.2. Plane wave method in photonic crystals study (PWM) 


Computational modeling of PBG devices has traditionally been approached using 
plane wave expansion techniques. We can take advantage of the unique properties 
of PCs with a calculation method necessary to determine how light will propagate 
through a particular periodic structure. For a periodic dielectric structure we can 
find eigenfrequencies (the allowed frequencies) for light propagation in all crystal 
directions being also able to calculate the field distributions in the crystal for any 
frequency of light. One of the most studied and reliable method is the PWM. It is 
simple enough to be easily implemented and it was used in some of the earliest 
studies of photonic crystals [10-12]. 


Modifying the conventional plane wave expansion method we can enlarge its 
application region making it suitable for solving the interface problem between 
vacuum and a semi-infinite PC system. This method has several advantages such 
as solving the problem of the ‘dielectric rods in air’ and ‘air rods in dielectric’, 
without any restriction on the shapes of the rods and position of cutting plane that 
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separating the semi-infinite PC region from the air region and we can use all 
information getting from the complete set of the eigenfunctions of the translation 
operator, including both the propagating and evanescent modes. So, we can 
analyze it easier and discuss the phenomena using the well-established knowledge 
of solid state physics. Additionally, we can easily isolate the finite size effects 
such as the resonance behavior of the transmittance curve caused by the finite 
thickness of the PC sample and accurately calculate the transmittance and 
reflectance for a very thick PC sample. 


This numerical method has produced successful results, but has disadvantages 
such as slow convergence and high demands for computing power and memory. 


An example of plane wave method implementation is given in the fourth part of 
this paper, where we considered two dielectric media with different dielectric 
constant and we represented the band diagram for a two-dimensional photonic 
crystal in a triangular lattice. The band structure was computed using PWM 
consisting in Fourier series expansion of the periodic functions and insertion of 
these series in the wave equation. By solving the problem we obtain a spectrum of 
eigenvalues (band diagram) and growth coefficients for Bloch eigenmodes [13]. 


2.3. Other methods 


2.3.1. Finite difference time domain method in photonic crystals study 
(FDTD) 


Another method developed and applied in the analysis and investigation of 
photonic crystals is the Finite Difference Time Domain method. It is a very 
general method which does not rely on any assumptions of propagating modes. 
FDTD method involves a direct discretization of Maxwell's equations, so that the 
evolution in space and time of a starting field (an incident wave) can be 
calculated. 


In [14] we can find a very helpful overview of the method with rich attention to 
implementation details and boundary conditions. A drawback of this method is the 
number of discretization points required for getting stability and accurate results 
from the algorithm. The time step should be less than the corresponding fraction 
of the smallest wave period in the problem and the spatial step size should be less 
than about a tenth of the smallest wavelength. 


This limits the size of the problems that can be handled on a given computer, the 
needed processing time being directly proportional to the number of grid points. 
We can calculate in this manner a band gap in a PC with an intuitively simple 
method. 
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We also can determine the transmission of waves in all directions through the 
crystal over the range of frequencies of interest. It involves calculations for many 
angles of incidence and many frequencies over a spatial region. 


In many practical application photonic crystals are considered materials with a high 
refractive index contrast (necessary in order to obtain an appreciable band gap) 
where waves having frequencies amid the band gap will penetrate for only a few 
periods, so we can determine its band. gap.only looking at the transmission through 
a relatively thin ‘block’ of PC material. 


Also, using a single calculation of FDTD algorithm we can obtain the results for all 
angles and frequencies. 


A drawback is that the method of obtaining frequency spectra by transforming a 
pulse response from the time domain calculations does not properly accounts for 
material dispersion which is specified in the frequency domain. 


However a method for incorporating material dispersion in FDTD calculations is 
known [14]. 

On the other hand, in some cases material dispersion may be neglected with respect 
to the large ‘geometric’ dispersion of the periodic PC itself [15]. 


We can determine the dispersion with FDTD method without performing a formal 
modal analysis, for a waveguide formed by removing a row of pillars from a 
photonic crystal consisting of 2D square array of square cross-section high-index 
pillars in air [4], [10-11], [16]. 


Also, we can determine the group velocity from the dependence of the frequency. 


The exited modes in the waveguide can be separated by spatial Fourier 
decomposition of the field pattern. 


The efficient coupling of waves into a waveguide in a photonic crystal is another 
practically important issue. Many phenomena can be modeled in PC structures with 
the FDTD method, some of them being the calculation of a photonic band gaps, 
resonator-based filters, diffraction effects and waveguide dispersion. 


There are limitation of the FDTD such as the difficulty of mcorporating material 
dispersion and the required computer resources. 


The illustration of a standard Cartesian Yee cell [17] used for FDTD can be seen in 
figure 2.1. 


We can visualize the electric field strength components (E,,E,, £. ) form the edges 
of the cube and the magnetic field strength components (H,,H,,H,) form the 
normals to the faces of the cube. 
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Fig. 2.1. The illustration of a standard Cartesian Yee cell. 


Using this model of discretization of Maxwell’s equations in time and space one 
can achieve the analyses over simple optoelectronic devices whose operation are 
based on the electromagnetic field propagation through periodic media. 


FDTD is an intuitive technique based on Maxwell’s equations, so one can easily 
understand how to use it or know what to expect from a fixed model. Being a time 
domain technique if the source we use is a broadband pulse then one can obtain a 
response over a wide range of frequencies. This fact can be used in applications 
when a broadband result is desired or when the resonant frequencies are not 
known. On the other hand, models with long features are difficult to model in 
FDTD, because of the large computational domain required. 


2.3.2. Multipole method in photonic crystals study (MM) 


Multipole methods have evolved as an important class of theoretical and 
computational techniques in the study of PCs. The MM is a numerical technique 
used as mode solver, a techniques that is particularly useful in finding modes of 
microstructured or photonic crystal fibers. It achieves high accuracy and rapid 
convergence with modest computational resources and it yields both the real and 
the imaginary parts of the mode propagation constant (provides information about 
losses), being used for full-vector modal calculations of microstructured fibers 
and other photonic crystal structures. 


The MM is implemented for circular inclusions in freely available software, but it 
can be extended to noncircular inclusions, too. It can deal with two types of 
microstructured optical fibers (MOF): air core MOF and a solid core MOF, 
surrounded by air holes. We can model systems with large numbers of inclusions 
and the computational overhead is further reduced for structures with discrete 
rotational symmetries by exploitation of the symmetry properties of the modes. 
Some convergence problems arise as the spacing between the inclusions decreases 
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because the method is limited to nonintersecting circular inclusions. Each 
dielectric boundary in the system is treated as a source of radiating fields in this 
method. The application of a certain field identity relates the regular field in the 
vicinity of any scatterer to fields radiated by other scatterers and external sources. 
The technique has evolved to become an important tool in the solution of dynamic 
problems involving both finite and periodic systems even if the origins of the 
method lie in the solution of electrostatic problems for periodic systems. 


Additionally, the MM can analytically preserve the symmetry and exploit the 
natural basis of functions for the scatterers, yielding important physical insight 
into the scattering processes. The propagation constant in the MM follows from 
the calculation when the frequency is used as an input parameter. If we are dealing 
with dispersive media this feature adds a significant advantage because the 
appropriate refractive indices are known from the outset. 


Essential features of MMs for periodic systems are lattice sums, which are sums 
of terms evaluated at each point of the lattice structure, taking into account its 
dependence on dimensionality and geometry. 


The MM results from considering the balance of incoming and outgoing fields and it 
allows solving the scattering problem associated with multipole inclusions [18-20]. 
We consider a single inclusion in the matrix (Fig. 2.1) centered at the origin of the 
coordinate system O. The field U(r,@ +27) =U(r,@)is periodic along the angular 
coordinate. We can expand U(r,@)in a Fourier series if we fix r: 
U(r, 0) = » f,(ne"” , where Fourier coefficients f(r) are regular functions of r. 


neZ 


Fig. 2.1. Single inclusion in the matrix where S, , are sources. 


If we consider an inclusion, the field reaching the corresponding region will be 
scattered. If we consider two inclusions, the incoming field for inclusion 1| results 
from the superposition of the field radiated from S, and the scattered field from 
inclusion 2. The same steps are followed for the incoming field for inclusion 2. 
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The MM is used in modeling of PCs and MOFs, to design complex/composite 
devices, such as interferometers, splitters, couplers and to study the radiation 
properties of sources embedded in such structures. Also, it is useful in studies of 
finite systems that involve the computation of the local density of states. 


2.3.3. Bloch wave method - Method of moments 


With MoM (Method of Moments) combined with Bloch wave expansion one can 
determine the band structure of photonic crystals. This approach is analogous to 
the spectral domain MoM method and is very efficient in terms of the number of 
plane wave needed for good convergence, being used for analyzing two- 
dimensional periodic structures. We can add that the field is expanded in terms of 
a plane wave spectrum [21] or in a set of eigen function modes (plane wave 
spectrum in 2D or a Bloch wave in 3D) and an integral equation is enforced on the 
surface of the scatterers in each unit cell. In the photonic crystal case the unit cell 
is three-dimensional and in the frequency selective surface case the unit cell is 
two-dimensional. 


The electric field E for ideal conducting structures admitting only electric current 
sources J is: 


B=—jky| A+2,v(v-A) (2.3.1) 


where A is the vector potential AA+k?A=-J and k is the free space 
propagation constant and it comes from the derivative in Maxwell’s equations. 


1 Sie Brt NG ancdi (Ca Brte) | _ £ 
ADE kx ty? — BES 2 - (2.3.2) 
i a,b,c a b Vc 


where we consider a cubic lattice in which a only depends on a, B only depends 
on b, y only depends on c and: 


a 


2 
a, =k, sing, cos dy = 


x 


B, =k, sin @, sing, ee ‘ 


y 


270C 
Y. =k, cos @y ae 


where k, =272//A, d,.. are the dimensions of the unit cell in the x, y, z 


directions, 4 is the effective wavelength in the crystal, g, and @, are the 


114 Dana Georgeta Popescu, Paul Sterian 


directions of propagation in spherical coordinates, k, represents the propagation 


constant in the periodic medium, /, represent unknown weighting coefficients 


WIGay.2= VE y,2)), 


1 a; a,b, Ail 
k? k? k? 
D 
Cx ie 1- 2 fs - is the Green function. 
aul Bro c 1- Ye 
2 k? 1” 


The (2.3.2) equation is very simple to implement and requires only the 3D Fourier 
Transform of the basis functions to be computed. Computing the band structure of 
a 3D photonic crystal with this method is as easy as computing transmission and 
reflection from a 2D periodic surface. Equation (2.3.2) can be used to compute 
bands in various types of doped and undoped photonic crystals [9-11]. 


We can assume values for (ky,@),@)) and then search for those values of k to 
compute bands of the crystal, which drive the determinant of the impedance 
matrix to zero. The (2.3.2) equation can be used to compute bands in various 
types of doped and undoped PCs [22-24]. 


3. Second harmonic generation (SHG) 


Solving nonlinear Helmholtz equations in combination with finite elements 
method, presented in section 2, and fixed point iterations reveals a nonlinear 
optical process called second harmonic generation (SHG). 


We can improve the SHG efficiency in nonlinear PCs by adjusting the 
fundamental wave frequencies and the second harmonic wave to the localized 
frequencies at the photonic band edge or by creating defects in the nonlinear 
materials of PCs. 


We consider a nonlinear material with periodic structure along the z direction, for 
0 < z< Dand the nonlinear structure is situated between two linear materials. The 
nonlinear Helmholtz equations [25-29] are: 


d°E, 2 2, (2) p* 
Fe + (kon) E, =-ky mE, E, (3.1) 
d°E 

; + (kon, eo = ~2ke yO EP (3.2) 


dz* 
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where E,, are the electric fields for the fundamental frequency x and for the 


second harmonic frequency 2x, E, is the complex conjugate, c is the speed of 
light in vacuum, k,) = @/c is the wave number in vacuum, n,, are the refractive 


indices of @ and 2 a, 7?) are two second order elements of the nonlinear 


susceptibility tensors for o and 2 o. 


The coupled nonlinear Helmholtz equation can be solved using the boundary 
conditions (for z = 0 and z = D). 


The structure we used is formed of 40 dielectric layers and the material is nonlinear. 
The refractive index alternates between value 1.000 and 1.428 for the fundamental 
wave and between 1.000 and 1.519 for the second harmonic wave and the layers are 
0.25 um and 0.35 um thick. The incident wave is the fundamental mode with 


E =1.0MV/m and the nonlinear susceptibility tensor vy” =0.1pm/V . In figure 
3.1a,b we represented the absolute value of the fundamental frequency field and 
second harmonic field for Q=0.592, where a huge amplification of the second 
harmonic field takes place. The conversion efficiency of the second harmonic 
should increase considerable because the fundamental frequency is increased with 
more than one order of magnitude in comparison with the peak value outside the 
structure. There are two maximum intensities of the SHW within each high index 
layer in contrast to the fundamental wave peak (the wavelength of the second 
harmonic signal is half of the pumped wave). 
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Fig. 3.1a. The absolute value of the fundamental frequency field (fundamental wave). 
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Fig. 3.1b. The absolute value of the second harmonic field (the second harmonic wave). 


In figure 3.2a,b are plotted the absolute values of the second harmonic field for a 
normalized frequency of 0.575 and 0.55, respectively. As we observed in figure 
3.1b, in picture plotted below are two maxima intensities of the second harmonic 
wave in each layer having 1.519 refractive index. 
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Fig. 3.2a. The absolute value of the second harmonic field (the second harmonic wave for Q=0.575). 
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We can see that the frequency decrease of the incident radiation results in 
significant changes of the shape and intensity of the resulting signal. As expected, 
for every incoming radiation field having different frequencies, the second 
harmonic is identified. Their spectral features closely rely on those of the 
incoming signals, supporting the accuracy of our calculations. 


x107(V/ son) 
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Fig. 3.2b. The absolute value of the second harmonic field 
(the second harmonic wave for Q=0.55). 


This periodic structure can be studied with FDTD method, too, a method 
described in the 2.3 section. The structure looks like the one in figure 3.3a. 
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Fig. 3.3a. One-dimensional photonic crystal lattice with 40 dielectric layers 
and four points of observation along the structure. 
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Figure 3.3a depicts the periodic structure we used to characterize the SHG in 
terms described earlier. Figure 3.3b illustrates the way the field propagates along 
the one-dimensional structure. 
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We can observe that the field is well collimated within y coordinates and it 
propagates freely. 
1. 


0.9 


0 10 30 40 


20 
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Fig. 3.3b. Field map and other results after simulation. 
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If we consider a photonic crystal with defect, a structure that we can see in figure 
3.4, changing the thickness from 0.304 um in 0.737 um of the 15" nonlinear layer 
(of 29 layers) we obtain a rise of the field in the defect. 


The linear layer is considered to be 0.351 um thick and the refractive index 1.000 
and 2.157 for the fundamental wave and 1.000 and 2.237 for the second harmonic 
wave. 


3,000 6,000 9,000 12,000 15,000 18,000 
a he tt 
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Fig. 3.4. One-dimensional photonic crystal lattice with 29 dielectric layers, three points of 
observation and one observation area along the structure. 


4. The plane wave method implementation in Matlab 


We consider two dielectric media with different dielectric constants | (air) and 
13 (GaAs). The lattice constant is 1.0. 


In figure 4.1 are represented calculated photonic band gaps for a two-dimensional 
photonic crystal consisting of cylinders with circular cross section and infinite 
length, in a triangular lattice. 


We consider plane propagation and two independent polarized states with E and H 
poles. We can observe a band gap situated between 0.216 um and 0.325 um for 
TE modes (continuous lines) and no band gap for TM modes (dashed lines). 
In conclusion, this structure has not a total band gap. 


To understand the plane wave method presented in section 2 we will expose the 
convergence problem and the dielectric function representation. 


To represent ¢,(7) in Fourier space we calculate Fourier coefficients of the 
inverse dielectric function 1/é,(7): 


nlG @)= S)ar eer i(G G).7) 


where G and G are two arbitrary vectors of reciprocal network, © is the 
primitive cell volume. 


Considering Maxwell’s equations we obtain for the magnetic field: 
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Vx and V-A(7)=0. 


2 
6,A)-(2) AG), vies 6, =0% 


E, (7 ) 


That is a standard eigenvalues problem, OH =AH. © y 1S an hermitic operator 
and has real eigenvalues. 


So, the magnetic field strength H becomes: 
H(F)=explik -7)H(*)%,, WF) = AF +R), A(F)= Si, exp (G-7) 
G 


where R is an arbitrary vector of the network, @, is the unit vector which is 
perpendicular over k and parallel with A. 
So, we obtain: 


H(7)= é, explik . FEW, elie 7)= SAG, Aexpli(k =) Pe ne 


G,,A 
where €,,,, is the unit vector perpendicular over LaG. nG,.A) is the 
amplitude component over the unit vector. 
Using Fourier expansion we obtain the matrix form of the PWM equation: 
A A ‘ “A 2) 

Yk +G|é+GJe(E-E 6° Se 76° |e _(@ Ie 

a —@6:@e eee |r c) | hye 
where e,- and e,. are polarization vectors. 


And for TM and TE modes the equation are: 
2) 
Yk +G|e+b(G-G)hg - (2) Ing 
Ye+che+a)"6-6),, (2), 


G 


To plot the TE and TM modes in the band diagram we first introduced in the 
program two dielectric media with different dielectric permittivity 1 and 13. 


We considered a two-dimensional periodic network with the network constant 
a=1um. In the elementary cell we have just one “atom” with the radius of 


R=0.48 a. Knowing the radius and the network constant, we calculated the 
filling factor. 
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The band diagram for a 2D PhC 


Frequency a/2mc 


0 02 04 06 08 1 12 14 
Wave vector ka/2n 


Fig. 4.1. The photonic band gaps for a two-dimensional photonic crystal. The red lines represents 
TM modes and the dashed ones TE modes. For TE modes we can observe a band gap for 0.216- 
0.325 frequency interval, and for TM modes we do not have a forbidden band. 


Next we introduced the direct lattice vectors 4d,,d,: d,=ay and 


d, =a(V3f+)/2 and the reciprocal lattice ~—_—_- vectors by, b,: 


— Bails ry ee 
Vaya, 2 V3a 


The direct lattice has a periodicity in the real space. The third part of the code asks 
us to introduce the n number. We considered a finite lattice with (2n+1)* 


elementary cells, even if the two-dimensional network is, theoretically, infinite. 
The “atoms” positions are given in the direct lattice by the position vectors R and 
in the reciprocal lattice by G. They are indexed by count variable and their 


number is equal with the number of the plane wave (2n +1), meaning that each 
atom is described by a plane wave. 


The bigger the number of the waves, the better the approximation. Knowing the 
dielectric permittivity lattice which depends of the position vectors of the 
reciprocal lattice G, we could indicate a k, trajectory on the edge of irreducible 
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Brillouin zone, passing through the high symmetry points I,X,M_ that 
correspond to k,,k,,k, vectors. Using the last equations for TE and TM modes 
we extend them to: 


TE: M=(real(k+G.')*real(k+G)+imag(k+G.')*imag(k+G)).*inv(eps2) 
TM: Me=abs(k+G.')*abs(k+G).*inv(eps2) 


The eig(M) function calculates.eigenvalues of M_ lattice, which are then sorted 
with sort(abs(eig(M))) function and at the end we chose only 10 eigenvalues with 
sqrt(abs(E(1:10))).. They are saved in variable freq(:,counter). 


We obtained ten energy bands because M lattice depends on the k wave vector. 
The band structure was finally obtained as depicted in Figure 4.1. 


5. Conclusions 


We have outlined some methods which may be used for characterizing PC. The 
techniques considered were Finite Element Method, which is a numerical 
technique used in finding approximate solutions of integral equations, or of PDE, 
Plane Wave Method, a method which consist in Fourier series development of the 
periodic functions and insertion of these series in the wave equation, Finite 
Difference Time Domain Method, an intuitive technique based on Maxwell’s 
equations which can be used in applications when we desire a broadband result or 
when we do not know the resonant frequencies, Multipole Method, which results 
from considering the balance of incoming and outgoing fields and with it one can 
solve the problem of scattering consisting of multipole inclusions and Bloch 
Wave Method which helps us determine the band structure of photonic crystals. 


Then the second harmonic problem was analyzed using a highly specialized 
program. We obtained the result that we expected, that meaning two maximum 
intensities of the second harmonic wave within each high index layer. At the end, 
we have plotted a band structure of a 2D PC using the Plane Wave Method. 
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